Geometric stabilization of extended S = 2 vortices in two-dimensional photonic 
lattices: theoretical analysis, numerical computation and experimental results 
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In this work, we focus our studies on the subject of nonlinear discrete self-trapping of S = 2 
(doubly- charged) vortices in two-dimensional photonic lattices, including theoretical analysis, nu- 
merical computation and experimental demonstration. We revisit earlier findings about S — 2 
vortices with a discrete model, and find that S = 2 vortices extended over eight lattice sites can 
indeed be stable (or only weakly unstable) under certain conditions, not only for the cubic nonlin- 

' earity previously used, but also for a saturable nonlinearity more relevant to our experiment with 

a biased photorefractive nonlinear crystal. We then use the discrete analysis as a guide towards 
numerically identifying stable (and unstable) vortex solutions in a more realistic continuum model 
with a periodic potential. Finally, we present our experimental observation of such geometrically ex- 
tended S — 2 vortex solitons in optically induced lattices under both self-focusing and self-defocusing 
nonlinearities, and show clearly that the S — 2 vortex singularities are preserved during nonlinear 

. propagation. 
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INTRODUCTION 



Over the last two decades, the study of Hamiltonian lattice systems, as well as of continuum models with periodic 
potentials has been a subject of increasing interest These systems arise in a diverse host of physical contexts, 
describing, e.g., the spatial dynamics of optical beams in coupled waveguide arrays or photorefractive crystals or 
optically induced photonic lattices in nonlinear optics Q , the temporal evolution of Bose-Einstein condensates (BECs) 
in optical lattices in soft-condensed matter physics [3], or the DNA double strand in biophysics among many others. 

One of the main thrusts of work in these directions has been centered around the investigation of existence and 
stability of localized solitary wave solutions. In two dimensions, such structures can be discrete solitons H, Q or 
discrete vortices (i.e., structures that have topological charge over a discrete contour) @, H[. Optically-induced 
photonic lattices in photorefractive crystals such as strontium barium niobate (SBN) have been used as an ideal 
platform for the observation of those predicted soliton structures. Indeed, the theoretical proposal [5| of such lattice 
solitons was followed quickly by their experimental realization in 2D induced lattices d, [To[ , subsequently leading to 



, the observation of a host of novel solitons in this setting, including dipole multipole [12], necklace [13[, and rotary 
[l4| solitons as well as discrete [l6] and gap [17] vortices. In addition to lattice solitons, photonic lattices have 
ON . enabled observations of other intriguing phenomena such as higher order Bloch modes [18], Zener tunneling [l9j], and 
localized modes in honeycomb [201 ], hexagonal [2l| and quasi-crystalline [22[ lattices, and Anderson localization [23[ 
^ , (see, e.g., the recent review [24] for additional examples). In parallel, experimental development in the area of BECs 
closely follows, with prominent recent results including the observation of bright, dark and gap solitons in quasi-one- 
dimensional settings [25], with the generation of similar structures in higher dimensions being experimentally feasible 
for BECs trapped in optical lattices 26|, [27| . 

Earlier experimental work on discrete vortex solitons [H, Gil nas mainly focused on vortices of unit topological 
charge (i.e., 5=1, with a 2tt phase shift circle around a discrete contour). However, more recently both in optics 
HQ and in BEC [13, EU (so far in the absence of the lattice for the latter case) , the study of higher charge vortices 
has been of interest. In particular, in the emerging area of hexagonal [2l|, [32[ and honeycomb [20[ lattices, it has 
been predicted [H, [34] and experimentally observed [35[ very recently that a higher-order vortex with topological 
charge 5 = 2 is more stable than a fundamental vortex with unit charge (5 = 1) when self-trapped with a focusing 
nonlinearity. 

The prototypical nonlinear dynamical lattice associated with the above systems is the so-called discrete nonlinear 
Schrodinger equation (DNLS) [36]. In the context of that model, it has been predicted that genuine 5 = 2 vortices 
turn out to be unstable in the case of a square lattice [7|, |8[]. This instability, however, can be avoided as proposed in 
two separate ways: a more intrusive one in which the central site of the contour is eliminated, introducing a defect 
therein [37[, as well as a less intrusive one involving the configuration of 5 = 2 vortex on a "rhomboidal" contour as 
8-site excitation in the square lattices, proposed solely on the basis of numerical observation as in [38[. Our aim in the 
present paper is to understand from a theoretical perspective the geometric stabilization (as we will call it) of 5 = 2 
vortices, and to illustrate its generic nature in discrete systems, by considering the case of optically induced photonic 



2 



lattices with a photorefractive nonlinearity (where there are some special intricacies of the theoretical analysis that we 
also illustrate). Based on our theoretical analysis of the discrete model, we then consider a more realistic continuum 
model of beam propagation with a periodic potential, illustrating that the S — 2 vortices geometrically extended to 
eight sites can indeed be linearly stable even in the continuum model. Finally, we corroborate these theoretical and 
numerical analyses with direct experimental observation of self-trapped S = 2 vortices whose topological charges are 
sustained during propagation throughout the nonlinear medium, which can be clearly contradistinguished from prior 
observations of breakup or charge-flipping of such high-order vortices [28] . 

The structure of our presentation is as follows. In section II we examine the model of the discrete nonlinear 
Schrodinger type with both Kerr and saturable nonlinearities, where our analytical findings are compared with nu- 
merical bifurcation results. In section III we consider the continuum model of beam propagation in photorefractive 
media with a periodic potential relevant to our experiment for the focusing and briefly also for the defocusing case. In 
section IV we present our experimental results. Section V concludes the paper with a number of interesting directions 
proposed for future studies. 

DISCRETE MODELS: ANALYSIS AND NUMERICS 

The DNLS with both Kerr and saturable nonlinearities can be written as: 

iu m ,n = -sA 2 U myri - A/"(|^ m ,n| 2 )^m,n (1) 

where Af(\u\ 2 ) = \u\ 2 for the typical DNLS with Kerr nonlinearity, while Af(\u\ 2 ) — —1/(1 + \u\ 2 ) in the saturable 
one (notice that both models share the same small amplitude limit). In the above, the overdot denotes the derivative 
with respect to the evolution variable, while A2'U m , n = Um+i,n + u m-i,n + ^™,n+i + Um,n-i — ^Um,n stands for the 
discrete Laplacian. We seek sationary solutions of the form ii m , n = exp(z/it)v m5n [notice that we denote the evolution 
variable as t, although in optics it represents the propagation distance z], obtaining the steady state equation: 

{/I, - eA 2 - N(\v m , n \ 2 )]v m , n = 0. (2) 

Our considerations will start at the so-called anti-continuum (AC) limit of e = 0, where it is straightforward to 
solve the steady state equations for a given \i. In the DNLS case, each site can be v m ,n = v^ ex P(^ m > n ) or v m,n = 
and in the saturable case, we have v m , n = — lexp(z# 

m,n) or ^m,n — 0- I n order to consider cases where the 
solutions have the same amplitude between our two examples, we will select, without loss of generality /i = 1 for the 
DNLS case, while \i = —1/2 for the saturable one. 

It is worthwhile to notice that in the AC limit, solutions can be obtained with arbitrary phase profiles. However, 
out of all the possible profiles, it is relevant to examine which ones may survive for e =^ 0. To do so, we expand the 
solution into a power series e.g., v = vo + ev\ + . . . , and write down the solution of the system order by order. To 
leading order the relevant equation for v\ will yield 



J(v ,e 




where J is the Jacobian of Eqs. ([2]) [with respect to the variables vq and Vq] evaluated at the solution for e = 0: 



7(v 2d(tfv)/dv -d{Nv)/dv* \ _ ( A 2 \ 

J[V ' £) ~ V -d{Nv)*/dv l-2d{Nv)*/dv* ) £ V A 2 / [6) 

For each of the excited sites (indicated by j) there corresponds a zero eigenvalue of the Jacobian with an eigenvector 
which vanishes away from the j-th block. Projecting to this eigenvector, as explained in 0, Q yields the Lyapunov- 
Schmidt conditions for the persistence of the solution. 

First, let us clarify the notation we will use in this paper. When the principal axes of a square lattice are normally 
oriented in horizontal and vertical directions, we name the structure a square vortex when it excites horizontal and 
vertical (nearest-neighbor) sites within a square contour, and likewise, a rhomboidal vortex when it excites diagonal 
(next-nearest-neighbor) sites in a rhomboidal contour (see Fig. [1]). We shall stick to such notation later for the 
continuum model even when the principal axes are diagonally oriented, as typically used in experiments with induced 
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FIG. 1: (Color online) The left panels show a typical example (for e — 0.1) of the geometrically stabilized rhomboidal S — 2 
vortex obtained from DNLS with Kerr (cubic) nonlinearity (top left is the real part, top right is its imaginary part, bottom 
left shows the amplitude and bottom right the phase distribution) . The right panel shows the dependence of the corresponding 
eigenvalues as a function of e. The solid lines show the full numerical results for the imaginary parts of the eigenvalues Xi 
as a function of the coupling constant e, while the two dashed lines (the lowest one of which is indistinguishable from the 
corresponding solid line) show the explicit analytical predictions in the text, which are in very good agreement with the 
corresponding numerics. 



lattices [15|, |16. 



For the case of the S — 2 square vortex, different conditions have been analyzed in 0, [33] 
and will not be presented here. Instead, we present here the case of the geometrically stabilized S = 2 rhomboidal 
vortex. In the latter case, if we denote the vortex by indices 1, ... ,8, starting from the top of the rhombus, then 
it is straightforward to see that the fundamental difference between the square and rhomboidal settings is that all 
interactions between adjacent sites in the rhombus are next nearest neighbor ones, and hence they come in not at the 
leading order in perturbative corrections (i.e., at 0(e)), but rather at 0(e 2 ), having the following form: 



g 2j +i = 2 sin(# 2 j+i - 2j ) + 2 sin(<9 2 j+i - #2.7+2) 

g 2j = 2 sin(<9 2 j - 2j -i) + 2 sm(6 2j - 2j+1 ) + sm(6 2j - 6 2j - 2 ) + sm(6 2j - 2j+2 ) 



(4) 
(5) 



where j = 1,2,3,4, with periodic boundary conditions (i.e., site corresponds to site 8, site 9 to site 1 etc.). Both 
in the DNLS, as well as in the saturable lattice case, the S = 2 solution with AO = tt/2 among adjacent sites (so 
that an accumulation of phase of 8 x tt/2 = 4tt occurs around the discrete contour for these extended 8-site vortices) 
naturally satisfies the above persistence conditions. 

However, the more crucial question is that of stability of the pertinent solutions. As initially illustrated in @,@], 
perhaps not surprisingly, the matrix that bears the relevant information is the Jacobian of the persistence conditions 
{M)j,k = dgj/dOk- To leading order, in the cubic case, the eigenvalues jj of this Jacobian are associated with the 
small eigenvalues of the full problem via 

X 2 j = 2e 7i . (6) 
The linear stability problem of solutions of Eq. (pQ) is then determined from the eigenvalue problem: 

J(v,e)r/> = i\o-i/>. (7) 

and a is 

/ 
-I 

Note that this system can equivalently be rearranged into 2-by-2 block formation. The discrete vortex is called 
spectrally unstable if there exist A and 1/? in the problem ([7j), such that Re(A) > 0. Otherwise, the discrete vortex is 
called spectrally stable. 

It is straightforward to observe that in the limit of s = 0, only the excited sites yield a set of N (i.e., as many as the 
number of such sites) pairs of null eigenvalues, while the non-excited sites yield eigenvalues with A = ±1 (which will 
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FIG. 2: (Color online) Plots of the imaginary part of the eigenvalues similarly to the right panel of Fig. [T] but obtained from 
DNLS with saturable nonlinearity. 



"become" the continuous spectrum of the problem); in the saturable case, this part of the spectrum is at ±1/2. Out 
of these N vanishing eigenvalue pairs in the AC limit, only one can be sustained at the origin for e ^ 0, due to the 
preserved phase invariance of the solution, while the remaining ones have to move off of the origin, as dictated by Eq. 
(|6]). Interestingly, in the case of the unstable square contour vortex, the phase change of tt/2 across sites renders zero 
all elements of the corresponding tridiagonal (in that case) Jacobian of the persistence conditions which, there, read 
gj = s'm(6j — Oj+i) ±sin(0j — Qj-i). This is the so-called "super-symmetric" case of [8] which needs to be treated at a 
higher order and, as a result, leads to higher order eigenvalues |8l. 1371] . These studies illustrated that in addition to a 
pair at and another pair of higher order, there is a quadruple pair A = ±^81, a single pair ±\/ a/80 ± Sei and a 



pair (responsible for the instability) A = ±y \/S0 — Se. The investigation of these works illustrated that this real pair 
was due to the interaction between the 4 sites adjacent to the central one of the vortex (to 2nd order, as mediated 



by the central site of the vortex). That is why the defect-induced stabilization worked in [37[ via its exclusion of this 
type of interactions. 

However, stabilization also ensues, in a geometric fashion, in the case of the rhomboidal 5 = 2 pattern. What 
happens in this case is that indeed the odd sites of the vortex still interact between them through the central site, 
however, now this interaction is geometrically "screened" by their lower order interaction with the even- numbered 
sites within the contour. As a result the instability is no longer mediated. More specifically, the analytical calculation 
in the case of the DNLS yields an 8 x 8 Jacobian, which for AO = tt/2 has 4 zero eigenvalues (which will become 
nonzero at a higher order), while the remaining 4 eigenvalues satisfy 

7i = -4sin 2 ^), (8) 

yielding, in addition to the phase- invariance induced persistent zero pair, 71,3 = —2 and 72 = —4. Because the 
interaction arises due to 2nd order neighbors, as shown in 0, H[, the relevant contribution is Xj = ±y^2jjS^ yielding 
a double pair Ai 5 3 = ±2ei and a single pair A2 = ±2y/2ei, in excellent agreement with the numerical results, as 
shown in Fig. [T] (especially, the ±2ei prediction can not be distinguished from its numerical counterpart). Notice 
that, additionally to the above 0(e) eigenvalue pairs, there are also five pairs of smaller eigenvalues, a double and two 
single ones, as well as one at 0, due to the persistence of the gauge symmetry. 

In the saturable case, we also observe numerically the same geometric stabilization effect, as is illustrated by the 
corresponding eigenvalues, for this case, of Fig. [2j However, here there is an interesting theoretical feature that is 
worth discussing. In the recent exposition of the saturable case in [39[, it was illustrated for nearest neighbor contours 
that the formula associating the eigenvalues of the full problem with those of the reduced Jacobian is: X 2 - = (l/2)7j£ 2 , 
due to the differences in the corresponding linearization operators. However, we observed numerically in the present 
(next-nearest-neighbor) setting that this formula no longer holds. In particular, it is found that the dashed lines in Fig. 
[2j which optimally match the numerical findings for small e, are given by A = ±2ei and A = ±\[2si^ namely they are 
less than their cubic counterparts by a factor of y/2 (rather than a factor of 2, as per the relation of [39[). Retracing 
step-by-step, the derivation of section 5 of [6], one can indeed theoretically identify this nontrivial difference, when 
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FIG. 3: (Color online) The complete bifurcation structure of the S = 2 square (thick, black) and rhomboidal (thin, red) 
vortices are presented in the top panels, with the power curve on the left and the stability curve on the right represented by 
the maximum real part of the linearization spectrum. The insets show close-ups near the saddle-node bifurcations that occur 
for both branches near the band edge. The bottom right panel is the lattice intensity pattern, where the eight vortex sites are 
marked by letters. The "square" vortex, as described in the discrete setting, is indicated by A-H, while the "rhomboidal" vortex 
is indicated by AB'CD'EF'GH'. The panels in the bottom left show the modulus, \U\ 2 , of prototypical solutions for the square, 
(6 S ), and rhomboidal, (b r - after destabilization close to the band edge), vortices as defined in the text. To the right of these 
are their respective linearization spectra and the insets in the top show the complex argument, or phase, while those to the 
bottom show the modulus in Fourier space (with axes in the standard horizontal and vertical orthogonal directions). 

repeating the relevant calculation for the saturable case. In particular, it is true (we have checked that this does hold 
true for saturable nonlinearities also in Id next-nearest-neighbor contours) that generally for saturable nonlinearities 
but next-nearest-neighbor contours, the small eigenvalues of the full problem are given by: 

A 2 c = e 2 M 2 c (9) 

where M. 2 denotes the next-nearest-neighbor Jacobian. In the cubic case, the corresponding formula (cf. (5.15) of 
0) has an extra factor of 2 in the right hand side. The resulting eigenvalues Xj = ±s^yj are compared to the full 
numerical results, yielding once again good agreement for small e in Fig. [2l 
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FIG. 4: (Color online) The modulus, \U\ 2 , of the remaining solutions marked in the bifurcation diagram in the top panel of Fig. 
[3] The linearization spectra are given in the insets in the top right. All solutions are for extended s=2 vortices with similar 
phase structures, but only the 8-site excitation in the form of square or rhomboidal vortices is stable. 



CONTINUUM MODELS: BIFURCATION ANALYSIS AND DYNAMICAL EVOLUTION 

For our consideration of the continuum problem, we use the non-dimensionalized version of the photorefractive 
model with saturable nonlinearity, as developed in detail in [l2|, [4o| , in the following form: 

iu = -A 2 u + AT(\u\ 2 )u (10) 

where A2 is now the 2D continuum Laplacian and Af(\u\ 2 ) — E/(l + I Q \ + \u\ 1 ). 

Here, u is the slowly varying amplitude of the probe beam normalized by the dark irradiance of the crystal I&, and 

'.. = W(2^)<^ (^), (ii) 

is a square optical lattice intensity function in units of Id. Here Iq is the lattice peak intensity, z is the propagation 
distance (in units of 2k\D 2 /7r 2 ), (x,y) are transverse distances (in units of D/tt), Eq is the applied DC field (in units 
of 7r 2 (/congZ) 2 r33) _1 ), D is the lattice spacing, fco = 27r/Ao is the wavenumber of the laser in the vacuum, Ao is the 
wavelength, n e is the unperturbed refractive index of the crystal for the extraordinarily polarized light, k\ = kon e , 
and r33 is the electro-optic coefficient for the extraordinary polarization. In line with the experiment, we choose the 
lattice intensity Iq = 5 (in units of Id). A plot of the lattice is shown in the bottom right panel of Fig. 02 also for 
illustrative purposes regarding the locations of the lattice sites excited by the vortex beam. Notice the lattice is now 
oriented diagonally, in accordance with the experiment, which implies that it is rotated by 7r/4 with respect to the x-y 
oriented lattice in the discrete model. As such, the rhombus in the present diagonal-oriented lattice is tantamount to 
the square in the x-y oriented lattice in discrete setting, and likewise the square would correspond to the rhombus. 
Nevertheless, we still call the vortex extended to cover ABCDEFGH sites a square vortex, and that extended to cover 
AB'CD'EF'GH' sites a rhomboidal vortex, according to the notation defined earlier. We choose the remaining physical 
parameters consistently with the experiment as 

D = 28 /im, A = 0.5 /im, n e = 2.3, r 33 = 280 pm/V. 

Thus, in the numerical results presented below, one x or y unit corresponds to 8.92 /im, one z unit corresponds to 
4.6 mm, and one E unit corresponds to 10.17 V/mm in physical units. We also set E = 13.76 which corresponds in 
dimensional units to 140 V/mm. 

It is well known that this model admits solutions of the form u(x,y,t) — U{x,y)e L ^ t for propagation constants 
in the semi-infinite gap, computed here as f3 > —5.632. We compute continuations of solutions in this region, as 
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FIG. 5: (Color online) The modulus, \u(t)\ 2 , of the dynamical evolution of slightly perturbed solutions from the bottom 
left of Fig. [3] are presented above by characteristic density isosurfaces of half-max the initial amplitude, i.e. D(x,y,t) = 
{(x,y,t);\u(x,y,t)\ 2 = (l/2)max x , y \u(x,y,0)\ 2 }. The initial condition is u(0) = U s (l + 0.0bUN[-l, 1]), where U s is the 
corresponding exact solution and UN[— 1, 1] is a uniform random variable in [—1, 1]. For the solutions of Fig. [3] the left image 
depicts the very mild instability of the unstable solution b s very close to the band edge, while the right image displays the 
strong instability of the configuration b r . 

shown in the top panels of Fig. [3] for both the square (thick, black) and geometrically stabilized rhomboidal (thin, 
red) S=2 configurations [43[. The rhomboidal vortex (with square contour in the present diagonal lattice) is clearly 
stable through most of the interval, in consonance with our calculations in the discrete problem. On the other 
hand, it is somewhat surprising to note that, for (3 sufficiently far from the linear spectrum, the square vortex (with 
rhomboidal contour in present diagonal lattice) can also be stable. There is clearly a bifurcation of the expected 
real eigenvalue pair through the origin, however, at (3 « —4.85. This confirms the prediction of the discrete model, 
although at the same time it indicates that its results should not be expected to be uniformly valid within the 
semi-infinite gap. The prototypical vortex solutions (b s ,b r ) are presented in the bottom left panels of Fig. El along 
with phases and corresponding Fourier space profiles (insets), as well as linearization spectra (right). The remaining 
solutions associated with the continuation of the relevant branches are presented together with linearization spectra 
insets in Fig. [U One can observe that similarly to what was shown for single-charge vortices in [41], the more 
strongly unstable counterparts of square and rhomboidal structures consist of patterns occupying additional wells of 
the periodic potential. 

The dynamics of perturbed solutions from the bottom left of Fig. [3] are presented in Fig. [5] by characteristic de nsity 
isosurfaces of half-max the initial amplitude, i.e. D(x,y,t) = {(#, y : t); \u(x, y, t)\ 2 = (l/2)max^ 52y |ii(x, y, 0)| 2 } [44j . 
The left image depicts the very mild instability of the unstable configuration denoted above as b s very close to the 
band edge, while the right image displays the strong instability through a real eigenvalue of the configuration denoted 
by b r . The mild instability of the former does not manifest itself until after t = 300, while the strong instability of 
the latter manifests itself by t = 100. Also, note that a four site breathing structure appears to persist in the latter 
case. 
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FIG. 6: (Color online) Density Profile (with phase inset) and linearization spectrum for the case of a defocusing nonlinearity 
for (3 = 9.2 (a,c) and 9.9 (b,d). The latter is very close to the band edge which in this case is at (3 « 10. 



Defocusing case 

For completeness, we consider the defocusing case briefly as well. The theoretical predictions for the discrete 
defocusing model are tantamount to the ones for the focusing case because the staggering transformation Vm^n = 
(— l) m+n i; mjn , which takes a solution v of the focusing problem to a solution v of the defocusing problem, does not 
affect next-nearest-neighbor configurations and takes the nearest-neighbor 8-site 5 = 2 vortex to an equivalent 5 = — 2 
(similarly to what happens for the nearest neighbor 4-site 5 = 1 vortex). 

In the continuum version of the model with the saturable nonlinearity, we use the transformation E = —E. The 
configurations now live shifted by one half period of the lattice to the right (i.e. each letter is shifted to the nearest 
minima to the right in the bottom right panel of Fig. [3]). The linear spectrum shifts and the localized solutions 
are now found within the first band-gap (as opposed to the semi-infinite gap for the focusing case). Similarly to the 
discrete model, again the principal predictions persist, as seen in Fig. [U I.e., there is a real pair of eigenvalues close 
to the band-edge for the square configuration and again there is a large stability region for the rhomboidal one. One 
difference is that the square configuration is now always unstable in the entire first bandgap due to complex quartets 
of eigenvalues (see Fig. [6fa) as an example). An example of this configuration after the real pair bifurcates through 
the origin is given in Fig. EJb). Rhomboidal solutions from before (c) and after (d) the destabilization close to the 
band edge are also shown in Fig. [6] along with their linearization spectra. These solutions again disappear near to 
the band edge (via saddle- node bifurcations) and do not bifurcate from linear modes. 



EXPERIMENTAL RESULTS 

In our experiment, we use a setup similar to that used in [l5[ for the observation of fundamental (5=1) discrete 
vortex solitons. The square lattice is induced in a biased photorefractive crystal (SBN:60 5 x 10 x 5 mm 3 ) by a 
spatially modulated partially coherent laser beam (A = 488 nm) sent through an amplitude mask. When the mask 
is appropriately imaged onto the input face of the crystal, the Talbot effect of the periodically modulated laser beam 
is suppressed by using a diffraction element, so the lattice intensity pattern remains invariant during the propagation 
throughout the crystal. The double-charged (S=2) vortex beam is generated by sending a coherent laser beam through 
a computer generated vortex hologram. In the experiment, the lattice beam is ordinarily-polarized and the vortex 
beam is extraordinarily-polarized. Thus the lattice beam will undergo nearly linear propagation while the vortex 
beam will experience a large nonlinearity due to the anisotropic photorefractive property. The input / output intensity 
patterns of the vortex beams are monitored with CCD cameras. In addition, the vortex beam exiting the crystal is 
also sent into a Mach-Zehnder interferometer for phase measurement, as needed. 

To observe self-trapping of doubly-charged vortex solitons, the donut-like vortex beam is expended and launched 




FIG. 7: (Color online) Experimental observation of a double-charge discrete vortex soliton extended to eight lattice sites under 
self-focusing nonlinearity. Top panels: (a) Input lattice beam pattern where the circle indicates the location of the S=2 vortex 
beam at input, (b, c) interference patterns of the input vortex with an inclined plane wave (b) and with a spherical wave (c), (d) 
output of the S=2 vortex after 10 mm of linear propagation through the crystal, showing the breakup into two well-separated 
S=l vortices. Bottom panels: (e) output vortex pattern at a low bias field, (f) self-trapped S=2 vortex pattern at a high 
bias field, (g) interference pattern of the vortex soliton with an inclined plane wave (zoomed), and (h) nonlinear output of the 
double-charge vortex without the lattice. Small circles in the inter ferograms mark the locations of vortex singularities. 
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FIG. 8: (Color online) Experimental observation of an extended S=2 gap vortex soliton in photonic lattices with self-defocusing 
nonlinearity. (a) Similar to Fig. [3(a) but the vortex location is different, (b) output vortex pattern at a low bias field (- 
700V/cm), (c) self-trapped vortex pattern at a high bias field (-1.5kV/cm), (d) zoom-in interferogram as in Fig. 7(g), and (e) 
the Fourier-space spectrum of the S=2 vortex soliton. 



into the lattice such that the vortex ring covers eight lattice sites [indicated by a blue circle in Fig. [71(a)], while 
the vortex core is overlapping with the central non-excited site. This arrangement corresponds to the square vortex 
configuration (ABCDEFG) illustrated in Fig. El The phase singularity of the input vortex is identified from two 
different interferograms shown in Figs. [3(b, c) (zoomed in so as to see the fringes more clearly). Two- fork fringes 
[Fig. [3(b)] and two-start spirals [Fig. 0(c)] in the central region of the interferograms clearly show the S=2 phase 
dislocations. As expected, the S = 2 vortex breaks up into two S=l vortices during linear propagation through the 
homogenous medium (i.e., without the waveguide lattice), as shown in Fig. [3(d). This is the natural outcome of 
the topological instability [42J. When a dc field is applied along the crystalline c-axis, the SBN crystal turns into a 
self- focusing medium 0, 0, Ho, Qj], EH EH . Under a proper strength of the nonlinearity, the vortex beam evolves into 
aS = 2 vortex soliton. Typical results are shown in the bottom panels of Fig. [3 for which the lattice period is about 
25/im and the vortex-to-lattice intensity ratio of the beams is about 1:4. At a low bias field of 800V/cm, the output 
vortex exhibits typical discrete diffraction covering many lattice sites [Fig. [3(e)]. However, when the bias field is 
increased to about 1.6kV/cm, the vortex beam self-traps into aS = 2 vortex soliton, covering mainly the eight sites 
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excited at the input [Fig. [7(f)]. In order to identify the phase structure of the nonlinear localized state, a tilted broad 
beam (akin to a quasi-plane wave) is sent to interfere with the vortex soliton at the output. It can be seen clearly [Fig. 
[7(g)] that the two forks remain in the center, and their bifurcation directions remain also unchanged as compared 
with the input vortex beam [Fig. [7(b)]. This provides a direct evidence for the formation of S — 2 high-order discrete 
vortex solitons with preserving phase singularities. We emphasize that such solitons are generated under the 8-site 
excitation, as predicted in theory, and they are quite different from the 4-site excitation which leads to dynamical 
charge flipping or disintegration of the topological charge, as observed in our previous experiment [28]. Without 
the lattice, the 5 = 2 vortex becomes unstable and breaks up into soliton filaments due to azimuthal modulation 
instability under the same bias conditions [Fig. [7(h)]. This indicates that, in an optically induced square lattice, the 
S=2 vortex can be geometrically stabilized and self-trapped, in good agreement with above theoretical and numerical 
results. 

Finally, Fig. [8]shows our experimental results on self-trapping of extended S=2 vortex in the photonic square lattices 
induced with a self-defocusing nonlinearity. Notice that the location of the vortex ring [illustrated by the blue circle in 
Fig. [8(a)] is different from that in Fig. [7(a) in order to have the extended 8-site excitation, now that the waveguides are 
located at the intensity minima (rather than maxima) of the lattice-inducing beam under self-defocusing nonlinearity 
[291 ]. In comparison with very recent study of S=2 discrete vortex solitons in hexagonal photonic lattices [33l. l34j . we 
found that the S=2 vortex can remain self-trapped [Fig. [8(c)] and maintain its singularities [Fig. [8(d)] also under the 
self-defocusing nonlinearity. In addition, differently from the more localized 4-site excitation of the S=2 vortex which 
evolves into a self-trapped quadrupole-like structure [291] , here the vortex phase structure remains. Furthermore, the 
Fourier-space spectrum [Fig. [8(e)] does not concentrate near the four M-points of the first Broullion zone, confirming 
the above numerical finding that the complex mode structures of such high-order vortex solitons do not bifurcate 
from the edge of the Bloch band. 

CONCLUSIONS AND FUTURE CHALLENGES 

In the present work, we have studied the geometric stabilization of S=2 vortices through the presence of next- 
nearest neighbor interactions in 2D square lattices. Different orientations for multi-site excitation of the vortices 
were studied and compared with previous work where the geometric stabilization is absent and nearest-neighbor 
interactions mediated through the central site leading to instability. Based on the results obtained from the standard 
cubic nonlinear Schrodinger lattice model, we extended our study to the case of the discrete model with a saturable 
nonlinearity, revealing some subtleties with respect to next-nearest neighbor eigenvalue calculations. Nevertheless, the 
principal analytical and numerical observations persisted therein. These findings were further studied by computations 
in the full continuum model, and we found that both the square and rhomboidal orientations of the extended vortex can 
lead to stable S=2 vortex solitons under appropriate conditions. Lastly, we corroborated these analytical and numerical 
results with experimental observations in optically induced photonic lattices. In a 10 mm long photorefractive crystal, 
we demonstrated that a double-charge vortex can maintain its singularity during nonlinear propagation in square 
lattices under both self-focusing and -defocusing nonlinearities. 

It would be interesting to experimentally study the dynamics of higher charge vortices, such as vortices of topological 
charge S = 3, for which the discrete model theory predicts potential stability, as well as to consider using crystals with 
longer propagation distances such that some of the above-predicted instability phenomena could be observable. On 
the other hand, from a theoretical perspective, a detailed understanding of vortices in more complex lattice settings 
such as superlattices and quasi-crystalline lattices could be important for relevant studies in other nonlinear systems 
involving vortices and periodic potentials. 

This work was supported by the 973 program, NSFC, PCSIRT, NSF, AFOSR and the Alexander von Humboldt 
Foundation. We thank C. Lou, L. Tang and J. Yang for discussions and assistance. 
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